clear all
close all
clc

global beta lambdaH lambdaEU A AgeDiffCoeff sigma OutsideCoeff delta gamma0 gamma1
global lambdaL lambdaEducM lambdaAgeM lambdaNation lambdaNatMen lambdaEducF lambdaNorth lambdaLIW lambdaLIW2 lambdaGIW lambdaGIW2 lambdaHNI
global AgeParam NationParam alpha Glob THETA Param SURPLUS VecTotSingles stderr
global indactive indA indbeta indlambda indlambdaH indsigma indlambdaH2 indOutside indAgeParam indNation inddelta indgamma
global dimS dimProv dimNum dimAge dimEduc dimNation dimHouse ageGrid educGrid typicalSurplus NorthSouth
global Match1 Match2 SinglesM1 SinglesF1 SinglesM2 SinglesF2 Weight ObsFertility
global PredictMarriage1 PredictMarriage2 PredictSinglesM1 PredictSinglesM2 PredictSinglesF1 PredictSinglesF2 ProvDistribMen ProvDistribWomen SimFertility
global CoeffDivorce gamma
global mincrit first 
global Italian ItalianLowLow ItalianLowHigh ItalianHighLow ItalianHighHigh ItalianMen ItalianWomen NonItalian Nationality IndMarriage

warning('off','all')

stderr=0;
mincrit=1e15;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Paramters defining the markets
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dimS=600;
dimNum=3;
ReducSize1=1;
ReducSize2=1;
LowParam=0;
dimSimul=2;
dimProv=21;
dimAge=6;
dimEduc=2;
dimNation=8;
dimHouse=2;
ageGrid=[22;27;32;37;42;50];
educGrid=[10;15];
nationGrid=kron(ones(dimAge*dimEduc,1),(1:dimNation)');
educGrid2=kron(educGrid,ones(dimAge*dimNation,1));
ageGrid2=kron(kron(ones(dimEduc,1),ageGrid),ones(dimNation,1));
Italian=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==1)&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])==1);
ItalianLowLow=Italian.*(repmat(educGrid2,[1 dimAge*dimEduc*dimNation])==10)&(repmat(educGrid2',[dimAge*dimEduc*dimNation 1])==10);
ItalianLowHigh=Italian.*(repmat(educGrid2,[1 dimAge*dimEduc*dimNation])==10)&(repmat(educGrid2',[dimAge*dimEduc*dimNation 1])==15);
ItalianHighLow=Italian.*(repmat(educGrid2,[1 dimAge*dimEduc*dimNation])==15)&(repmat(educGrid2',[dimAge*dimEduc*dimNation 1])==10);
ItalianHighHigh=Italian.*(repmat(educGrid2,[1 dimAge*dimEduc*dimNation])==15)&(repmat(educGrid2',[dimAge*dimEduc*dimNation 1])==15);
ItalianMen=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==1)&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])~=1);
ItalianWomen=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])~=1)&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])==1);
NonItalian=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])~=1)&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])~=1);
ForeignWomen=repmat(nationGrid',[dimAge*dimEduc*dimNation 1])~=1;
EUNEW=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==3)|(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])==3);
Italian=repmat(Italian,[1 1 dimProv]);
ItalianLowLow=repmat(ItalianLowLow,[1 1 dimProv]);
ItalianLowHigh=repmat(ItalianLowHigh,[1 1 dimProv]);
ItalianHighLow=repmat(ItalianHighLow,[1 1 dimProv]);
ItalianHighHigh=repmat(ItalianHighHigh,[1 1 dimProv]);
EUNEW=repmat(EUNEW,[1 1 dimProv]);

IndMarriage=zeros(96,96,dimNation,dimNation);

for i=1:dimNation
    for j=1:dimNation
        IndMarriage(:,:,i,j)=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==i)&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])==j);
    end
end
IndMarriage=logical(IndMarriage);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Position of parameters in vector THETA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
indbeta=1:6;
indsigma=7;
indlambda=8:20;
indlambdaH=21:23;
indA=24:27;
indlambdaH2=28:30;
indOutside=31:53;
indAgeParam=54:88;
indNation=89:123;
inddelta=124:127;
indgamma=128:129;
indactive=indgamma;
DefaultCoeff=0.25;


% Weight for the graphs
AgeMen=repmat(ageGrid2,[1 dimAge*dimEduc*dimNation dimProv]);
AgeMen=(AgeMen-20)*6;
AgeWomen=repmat(ageGrid2',[dimAge*dimEduc*dimNation dimProv]);
AgeWomen=(AgeWomen-20)/5;
first=1;

% Provinces:
% (7) Ancona, (8) Bari, (9) Bergamo, (10) Bologna, (11) Brescia, (12) Cagliari   
% (13) Catania, (14) Catanzaro, (15) Firenze,(16) Genova, (17) Laquila, (18) Milano,     
% (19) Napoli, (20) Padova, (21) Palermo, (22) Perugia, (23) Potenza, (24) Roma       
% (25) Salerno, (26) Torino, (27) Trento, (28) Trieste, (29) Venezia    
ProvinceIndexMarriage=[7:29];
ProvinceIndexSingle=ProvinceIndexMarriage-3;


%% initialisation of parameters
seed=211187;
alpha=0.965;   % fraction of italians that do not meet foreigners. Roughly 1/2 Italians, 1/2 foreigners
% alpha=0.97;
A=[1;2;0;0];
AgeParam=ones(dimAge*dimAge-1,1);
NatioParam=zeros(sum(sum(triu(ones(dimNation,dimNation),1))),1);

beta=[unconstrain(0.5,0,1);unconstrain(0.7,0,1)];     % penalty for not having access to legal work, North, South
lambdaNatMen=0;        % Nationality and gender
lambdaEducF=0;        % Gender and education
lambdaAgeM=0;         % Nationality and age for men 
lambdaL2=0;        % genetic distance
lambdaG=0;        % genetic distance
lambdaG2=0;        % genetic distance
lambdaA=0;        % linguistic distance for Italian wives
lambdaEU=0.05;     % cultural distance of EU citizens to Italians
lambdaMW=0.2;     % european men and women from poorer countries
sigma=0.2;          % variance of the love component
lambdaH=zeros(3,1);  % effect of housing on surplus
AgeDiffCoeff=[4;30;17];
AgeDiffCoeff=[4;2;0.1];
OutsideCoeff=zeros(25,1);
Policy=0;
Period=1;
Num=1;
CoeffDivorce=0;
%% Distribution of characteristics by province: MARRIAGES
% order is education, age, origin
% Columns are:
% (1) Male education, (2) Male age, (3) Male origin, (4) Female educ, 
% (5) Female educ, (6) Female origin
% Provinces:
% (7) Ancona, (8) Bari, (9) Bergamo, (10) Bologna, (11) Brescia, (12) Cagliari   
% (13) Catania, (14) Catanzaro, (15) Firenze,(16) Genova, (17) Laquila, (18) Milano,     
% (19) Napoli, (20) Padova, (21) Palermo, (22) Perugia, (23) Potenza, (24) Roma       
% (25) Salerno, (26) Torino, (27) Trento, (28) Trieste, (29) Venezia    
% Missing here Campobasso, Verona

dimProv=length(ProvinceIndexMarriage);

marriage_before=csvread('All_Marriages_Before.csv');
marriage_before=marriage_before(:,ProvinceIndexMarriage);
Match1=reshape(marriage_before,dimEduc*dimAge*dimNation,dimEduc*dimAge*dimNation,dimProv);
% Rows are men, columns women
Match1=permute(Match1,[2 1 3]);
marriage_after=csvread('All_Marriages_After.csv');
marriage_after=marriage_after(:,ProvinceIndexMarriage);
Match2=reshape(marriage_after,dimEduc*dimAge*dimNation,dimEduc*dimAge*dimNation,dimProv);
Match2=permute(Match2,[2 1 3]);
%% Provinces kept in analysis
%  (18) Milano, (26) Torino, (10) Bologna, (16) Genova, (15) Firenze, (22) Perugia, 
% (24) Roma, (19) Napoli, (8) Bari (21) Palermo (7) Ancona (12) Cagliari (14) Catanzaro
% (17) Laquila (23) Potenza (27) Trento (28) Trieste (29) Venezia (11)
% Brescia  (20) Padova   (25) Salerno
% Missing:  Aosta Campobasso
indMarket=[18 26 10 16 15 22 24 19 8 21 7 12 14 17 23 27 28 29 11 20 25]-6;

Match1=Match1(:,:,indMarket);
Match2=Match2(:,:,indMarket);
correct_before=[1 1;2 1;3 1.52;4 1.78;5 1.74;6 1.63;7 1.12;8 1];  % Comes from missing_marriages2001.dta
correct_after=[1 1;2 1;3 1.36;4 1.55;5 1;6 1.22;7 1;8 1];  % Comes from missing_marriages2011.dta

mat_before=ones(dimEduc*dimAge*dimNation,dimEduc*dimAge*dimNation,length(indMarket));
mat_after=ones(dimEduc*dimAge*dimNation,dimEduc*dimAge*dimNation,length(indMarket));
for i=1:dimNation
    ind=repmat((nationGrid==i)&(nationGrid'==i),[1 1 length(indMarket)]);
     mat_before(ind)=correct_before(i,2);
     mat_after(ind)=correct_after(i,2);
end

Match1=round(Match1.*mat_before,0);
Match2=round(Match2.*mat_after,0);


%% SINGLES
% (1) education, (2) age, (3) origin 
% Provinces:
% (4)  Ancona ,(5)  Bari, (6)  Bergamo,(7)  Bologna,(8)  Brescia       
% (9)  Cagliari,(10) Campobasso,(11) Catania,(12) Catanzaro,(13) Firenze       
% (14) Genova,(15) LAquila,(16) Milano,(17) Napoli,(18) Padova,(19) Palermo       
% (20) Perugia,(21) Potenza,(22) Roma,(23) Salerno,(24) Torino,(25) Trento        
% (26) Trieste,(27) Venezia,(28) Verona        
Singles_Men_Before=csvread('Singles_Men_Before.csv');
Singles_Men_After=csvread('Singles_Men_After.csv');
Singles_Women_Before=csvread('Singles_Women_Before.csv');
Singles_Women_After=csvread('Singles_Women_After.csv');
ProvinceIndexSingle=4:28;
Singles_Men_Before=Singles_Men_Before(:,ProvinceIndexSingle);
Singles_Men_After=Singles_Men_After(:,ProvinceIndexSingle);
Singles_Women_Before=Singles_Women_Before(:,ProvinceIndexSingle);
Singles_Women_After=Singles_Women_After(:,ProvinceIndexSingle);

% (16) Milano, (24) Torino, (7) Bologna, (14) Genova, (13) Firenze, (20) Perugia, 
% (22) Roma, (17) Napoli, (5) Bari (19) Palermo
% (4) Ancona (9) Cagliari (12) Catanzaro
% (15) Laquila (21) Potenza (25) Trento (26) Trieste (27) Venezia (8)
% Brescia  (18) Padova   (23) Salerno
indMarket=[16 24 7 14 13 20 22 17 5 19 4 9 12 15 21 25 26 27 8 18 23]-3;
%indMarket=indMarket(1:4);
Singles_Men_Before=Singles_Men_Before(:,indMarket);
Singles_Men_After=Singles_Men_After(:,indMarket);
Singles_Women_Before=Singles_Women_Before(:,indMarket);
Singles_Women_After=Singles_Women_After(:,indMarket);
dimProv=length(indMarket);

ProvDistribMen_before=cumsum(Singles_Men_Before)./sum(Singles_Men_Before);
ProvDistribWomen_before=cumsum(Singles_Women_Before)./sum(Singles_Women_Before);
RatioSingleMF_before=sum(Singles_Men_Before)./sum(Singles_Women_Before);
ProvDistribMen_after=cumsum(Singles_Men_After)./sum(Singles_Men_After);
ProvDistribWomen_after=cumsum(Singles_Women_After)./sum(Singles_Women_After);
RatioSingleMF_after=sum(Singles_Men_After)./sum(Singles_Women_After);

SinglesM1=Singles_Men_Before+sum(permute(Match1,[1 3 2]),3); % dEd*dAge*dO dProv
TotalM1=sum(SinglesM1);
SinglesM1=Singles_Men_Before./TotalM1;
SinglesF1=Singles_Women_Before+sum(permute(Match1,[2 3 1]),3); % dEd*dAge*dO dProv
TotalF1=sum(SinglesF1);

Weight=TotalM1+TotalF1;
Weight=Weight/sum(Weight);

SinglesF1=Singles_Women_Before./TotalF1;
SinglesM2=Singles_Men_After+sum(permute(Match2,[1 3 2]),3); % dEd*dAge*dO dProv
TotalM2=sum(SinglesM2);
SinglesM2=Singles_Men_After./TotalM2;
SinglesF2=Singles_Women_After+sum(permute(Match2,[2 3 1]),3); % dEd*dAge*dO dProv
TotalF2=sum(SinglesF2);
SinglesF2=Singles_Women_After./TotalF2;

auxTotalF1=repmat(TotalF1',[1 dimEduc*dimAge*dimNation dimEduc*dimAge*dimNation]);
auxTotalF1=permute(auxTotalF1,[2 3 1]);
auxTotalF2=repmat(TotalF2',[1 dimEduc*dimAge*dimNation dimEduc*dimAge*dimNation]);
auxTotalF2=permute(auxTotalF2,[2 3 1]);

Match1=Match1./auxTotalF1;
Match2=Match2./auxTotalF2;

%% Shares of singles by nationality in province
VecTotSingles=zeros(dimNation,length(indMarket),2);
auxBefore=Singles_Women_Before+Singles_Men_Before;
auxAfter=Singles_Women_After+Singles_Men_After;

for i=1:dimNation
    ind=repmat(nationGrid==i,[1 length(indMarket)]);
    VecTotSingles(i,:,1)=sum(auxBefore.*ind);
    VecTotSingles(i,:,2)=sum(auxAfter.*ind);
end
VecTotSingles(:,:,1)=VecTotSingles(:,:,1)./sum(VecTotSingles(:,:,1));
VecTotSingles(:,:,2)=VecTotSingles(:,:,2)./sum(VecTotSingles(:,:,2));

%% Divorces 
auxDivorceMoments=xlsread('moments_separations_3cohorts');
DivorceMoments=zeros(7,2,3);
DivorceMoments(:,:,1)=auxDivorceMoments(4:10,1:2);
DivorceMoments(:,:,2)=auxDivorceMoments(19:25,1:2);
DivorceMoments(:,:,3)=auxDivorceMoments(34:40,1:2);


%% Shadow economy
% from file shadow_economy.dta 2002 values
shadow=[8.1 9.5 8.5 12.2 9.5 13 13.1 22.2 18.2 21.9 10.5 17.2 26 13.6 19.3 8.6 10.7 8.8 8.1 8.8 22.2]/100;

% 1: Italian, 2: EU before 2004, 3: EU12 (new EU countries), 
% 4: Other_EU (European countries but not part of the EU), 5: Africa,
% 6: Asia,7: South_America, 8: OECD outside of Europe
% data=csvread('CulturalMeasures.csv',1,0);
% GeneticDist=reshape(data(:,3),8,8);
% GeneticDist=GeneticDist/max(GeneticDist(:));
% LinguisDist=reshape(data(:,4),8,8);
% LinguisDist=LinguisDist-min(LinguisDist(:));
% LinguisDist=LinguisDist/max(LinguisDist(:));
% ReligionDist=reshape(data(:,4),8,8);
% TotalDist=reshape(data(:,4),8,8);




data=csvread('homeownership_25provcsv.csv');
ProbHouse=data(:,7);
ProbHouse=reshape(ProbHouse,96,4);  % first two columns are men / women in first period, then men / women in second period.

%  Standard errors for marriage rates
auxsd=xlsread('moments_sd_marriages2'); % coefficients are starting on line 4. Col 1 are coeff, Col2 are st err of estimates
auxsd=auxsd(4:end,1); 
% Province effects
sdprov=[0;auxsd(26:47)]; 
sdprov=sdprov(1:21);
sdprov=repmat(sdprov,[1 dimEduc*dimAge*dimNation dimEduc*dimAge*dimNation]); % dProv (dAge*dNat*dEd) (dAge*dNat*dEd)
sdprov=permute(sdprov,[2 3 1]);

% Nationality effects
nataux=auxsd(1:23);
nateffect=zeros(dimNation,dimNation);
nateffect(logical(eye(dimNation)))=[nataux(1);nataux(16:22)];
nateffect(1,2:dimNation)=nataux(2:8);
nateffect(2:dimNation,1)=nataux(9:15);
nateffect=kron(ones(dimEduc*dimAge,dimEduc*dimAge),nateffect);

% Age effect
aux=kron(ones(dimEduc,1),kron((1:dimAge)',ones(dimEduc*dimAge*dimNation*dimNation,1)));
aux=reshape(aux,dimEduc*dimAge*dimNation,dimEduc*dimAge*dimNation)';
ageeffect=(aux==aux')*auxsd(24); 

% Educ effect
aux=kron((1:dimEduc)',ones(dimEduc*dimAge*dimNation*dimAge*dimNation,1));
aux=reshape(aux,dimEduc*dimAge*dimNation,dimEduc*dimAge*dimNation)';
educeffect=(aux==aux')*auxsd(25); 

SD_Match1=repmat(auxsd(49)+nateffect+ageeffect+educeffect,[1 1 dimProv])+sdprov;
SD_Match2=SD_Match1;
SD_Match1=SD_Match1+auxsd(48);
SD_Match1=exp(SD_Match1);
SD_Match2=exp(SD_Match2);

%% Fertility (DiD) Table 12 in fertility_moments.xlsx
ObsFertility=[0.527 0.527 0.105 0.0471;0.000463 0.000463 0.00735 0.0226];
ObsFertility(2,:)=ObsFertility(2,:).*[50 50 1 1];


%% Geographical positions of provinces
% 1 is North, 2 is South
NorthSouth=ones(dimProv,1);
NorthSouth(1)=1;

Glob.dimS=dimS;
Glob.dimSimul=dimSimul;
Glob.dimAge=dimAge;
Glob.dimEduc=dimEduc;
Glob.dimNation=dimNation;
Glob.dimProv=dimProv;
Glob.dimHouse=dimHouse;
Glob.ReducSize1=ReducSize1;
Glob.ReducSize2=ReducSize2;
Glob.LowParam=LowParam;
Glob.ageGrid=ageGrid;
Glob.educGrid=educGrid;
Glob.DivorceMoments=DivorceMoments;
% Glob.GeneticDist=GeneticDist;
% Glob.LinguisDist=LinguisDist;
Glob.Singles_Men_Before=Singles_Men_Before;
Glob.Singles_Women_Before=Singles_Women_Before;
Glob.Singles_Men_After=Singles_Men_After;
Glob.Singles_Women_After=Singles_Women_After;
Glob.shadow=shadow;
Glob.RatioSingleMF_before=RatioSingleMF_before;
Glob.RatioSingleMF_after=RatioSingleMF_after;
Glob.ProbHouse=ProbHouse;
Glob.seed=seed;
Glob.NorthSouth=NorthSouth;
Glob.Policy=0;
Glob.PolicyParam=0;
Glob.NationParam=NationParam;
Glob.DefaultCoeff=DefaultCoeff;
Glob.VecTotSingles=VecTotSingles;
iprov=1;

THETA=[beta;sigma;lambdaNatMen;lambdaAgeM;lambdaEducM;lambdaG2;lambdaL;lambdaG;lambdaL2;lambdaG2;lambdaL;lambdaG;lambdaL2;lambdaG2;lambdaEU;A(:);AgeDiffCoeff;OutsideCoeff;AgeParam;NationParam];


load THETA_20191030
load Sstore2
%%%%%%%%%%%%%%%%%%%%
THETA(indA(4))=0;
%%%%%%%%%%%%%%%%%%%%
theta=THETA;
if stderr==0
tic
crit=estimMarriage20191030(THETA(indactive));
toc
disp(['Criteria : ' num2str(crit)])

c1=corr([Match1(:) PredictMarriage1(:)]);
c2=corr([Match2(:) PredictMarriage2(:)]);
disp(['correlation Period 1: ' num2str(round(c1(1,2),2),2)])
disp(['correlation Period 2: ' num2str(round(c2(1,2),2),2)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TABLE 8
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delete('table8.tex')
diary('table8.tex')
disp('\multicolumn{7}{l}{A. Fertility}')
disp('\hline')
disp(['Native male - foreign female, DiD & ' num2str(ObsFertility(1,3)) '(' num2str(ObsFertility(2,3),3) ') &' num2str(SimFertility(1,3))  ' \\'])
disp(['Native male - foreign female, level & ' num2str(ObsFertility(1,1)) '(' num2str(ObsFertility(2,1),3) ') &' num2str(SimFertility(1,1))  ' \\'])
disp(['Native female - foreign male, DiD & ' num2str(ObsFertility(1,4)) '(' num2str(ObsFertility(2,4),3) ') &' num2str(SimFertility(1,4))  ' \\'])
disp(['Native female - foreign male, level & ' num2str(ObsFertility(1,2)) '(' num2str(ObsFertility(2,2),3) ') &' num2str(SimFertility(1,2))  ' \\'])
disp('\hline')
mat1=[squeeze(DivorceMoments(:,1,1)) squeeze(DivorceMoments(:,2,1)) squeeze(CoeffDivorce(:,1,1))];
mat2=[squeeze(DivorceMoments(:,1,2)) squeeze(DivorceMoments(:,2,2)) squeeze(CoeffDivorce(:,1,2))];
mat3=[squeeze(DivorceMoments(:,1,3)) squeeze(DivorceMoments(:,2,3)) squeeze(CoeffDivorce(:,1,3))];
str1=[repmat('&',[7 1]) num2str(round(mat1(:,1),4)) repmat('&(',[7 1]) num2str(round(mat1(:,2),4)) repmat(')&',[7 1]) num2str(round(mat1(:,3),4)) repmat('&',[7 1])];
str2=[num2str(round(mat2(:,1),4)) repmat('&(',[7 1]) num2str(round(mat2(:,2),4)) repmat(')&',[7 1]) num2str(round(mat2(:,3),4)) repmat('&',[7 1])];
str3=[num2str(round(mat3(:,1),4)) repmat('&(',[7 1]) num2str(round(mat3(:,2),4)) repmat(')&',[7 1]) num2str(round(mat3(:,3),4)) repmat('\\',[7 1])];
disp('\multicolumn{7}{l}{B. Separation}')
disp(['&\multicolumn{2}{c}{Before enlargements} & \multicolumn{2}{c}{During enlargements}& \multicolumn{2}{c}{After enlargements}\\'])
disp(['&Observed & Predicted&Observed & Predicted&Observed & Predicted\\'])
disp('\hline')
disp([str1 str2 str3])
disp('\hline')
diary off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TABLE 9
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delete('table9.tex')
diary('table9.tex')
disp('\begin{tabular*}{0.99\textwidth}{@{\extracolsep{\fill} } llll}')
disp('\hline')
disp('\hline')
disp([ 'No legal access ($\underline{\beta}$), No shadow economy & ' num2str(beta(1),3) '&(' num2str(round(setheta(indbeta(1)),3),3) ') \\'])
ind=[find(INDACTIVE==indbeta(1));find(INDACTIVE==indbeta(2))];aux=Sigma(ind,ind);aux=sqrt(sum(aux(:)));
disp([ 'No legal access ($\underline{\beta}$), shadow economy 100\% & ' num2str(beta(1)+beta(2),3) '&(' num2str(round(aux,3),3) ') \\'])
disp([ 'Legal access, native husband - foreign wife ($\bar{\beta}^{\protect\male}$), No shadow economy& ' num2str(beta(3),3) '&(' num2str(round(setheta(indbeta(3)),3),3) ') \\'])
disp([ 'Legal access, native wife - foreign husband ($\bar{\beta}^{\protect\female})$, No shadow economy& ' num2str(beta(5),3) '&(' num2str(round(setheta(indbeta(5)),3),3) ') \\'])
ind=[find(INDACTIVE==indbeta(3));find(INDACTIVE==indbeta(4))];aux=Sigma(ind,ind);aux=sqrt(sum(aux(:)));
disp([ 'Legal access, native husband - foreign wife ($\bar{\beta}^{\protect\male})$, shadow economy 100\%& ' num2str(beta(3)+beta(4),3) '&(' num2str(round(aux,3),3) ') \\'])
ind=[find(INDACTIVE==indbeta(5));find(INDACTIVE==indbeta(6))];aux=Sigma(ind,ind);aux=sqrt(sum(aux(:)));
disp([ 'Legal access, native wife - foreign husband ($\bar{\beta}^{\protect\female})$, shadow economy 100\%& ' num2str(beta(5)+beta(6),3) '&(' num2str(round(aux,3),3) ') \\'])
disp('\hline')
ind=[find(INDACTIVE==indOutside(1)):find(INDACTIVE==indOutside(7))];aux=Sigma(ind,ind);
D=[1+sum(THETA(indOutside(2:7)));repmat(THETA(indOutside(1)),[6 1])]/7;aux=sqrt(D'*aux*D);
disp([ 'Average outside value, before enlargement, North & ' num2str(mean(OutsideCoeff(1:7,1)),3) ' & (' num2str(round(aux,3),3)  ') \\'])
ind=[find(INDACTIVE==indOutside(1)) find(INDACTIVE==indOutside(8)):find(INDACTIVE==indOutside(21))];aux=Sigma(ind,ind);
D=[1+sum(THETA(indOutside(8:21)));repmat(THETA(indOutside(1)),[14 1])]/15;aux=sqrt(D'*aux*D);
disp([ 'Outside value, before enlargement, South & ' num2str(mean(OutsideCoeff(8:dimProv,1)),3) ' & (' num2str(round(aux,3),3)  ') \\'])
ind=[find(INDACTIVE==indOutside(1)):find(INDACTIVE==indOutside(7)) find(INDACTIVE==indOutside(22))];aux=Sigma(ind,ind);
D=[THETA(indOutside(22))+THETA(indOutside(22))*sum(THETA(indOutside(2:7))); ...
    THETA(indOutside(1))*1+sum(THETA(indOutside(2:7)));...
THETA(indOutside(22))*repmat(THETA(indOutside(1)),[6 1])]/7;aux=sqrt(D'*aux*D);
disp([ 'Outside value, after enlargement, North & ' num2str(mean(OutsideCoeff(1:7,2)),3) ' & (' num2str(round(aux,3),3)  ') \\'])
ind=[find(INDACTIVE==indOutside(1)) find(INDACTIVE==indOutside(22)) find(INDACTIVE==indOutside(8)):find(INDACTIVE==indOutside(22))];aux=Sigma(ind,ind);
D=[THETA(indOutside(22))+THETA(indOutside(22))*sum(THETA(indOutside(8:21))); ...
    THETA(indOutside(1))*1+sum(THETA(indOutside(8:21)));...
THETA(indOutside(22))*repmat(THETA(indOutside(1)),[15 1])]/15;aux=sqrt(D'*aux*D);
disp([ 'Outside value, after enlargement, South & ' num2str(mean(OutsideCoeff(8:dimProv,2)),3) ' & (' num2str(round(aux,3),3)  ') \\'])
% disp([ 'Outside value, after & ' num2str(mean(OutsideCoeff(1:dimProv,2)),3) ' \\'])
disp('\hline')
disp('&&Italian& EU15 & New & Other  & Africa & Asia &South & OECD\\')
disp('&& &  & EU&  Europe &  &  & America & \\')
Name1=['   ';'   ';'   ';'Men';'   ';'   ';'   ';'   '];
Name1='\parbox[t]{2mm}{\multirow{3}{*}{\rotatebox[origin=c]{90}{Men}}}';
Name1=repmat(Name1,[16 1]);Name1([1:7 9:16],:)=' ';

Name2=['Italian    ';'           ';'EU         ';'           ';'New EU     ';'           ';'Oth. Europe';'           ';'Africa     ';'           ';'Asia       ';'           ';'S America  ';'           ';'OECD       ';'           '];

N=round(Nationality,2);
ind=[find(INDACTIVE==indNation(1)):find(INDACTIVE==indNation(28))];aux=Sigma(ind,ind);
auxA=1/4/atan(1);
D=auxA./(1+THETA(indNation(1:28)).^2);
aux=sqrt(D.*D.*diag(aux));
Natse=zeros(dimNation,dimNation);
Natse(logical(triu(ones(dimNation,dimNation),1)))=aux;
Natse=Natse'+Natse;
D=[lambdaNatMen*log(Nationality(1,3:end).^lambdaNatMen); ...
    lambdaNatMen*Nationality(1,3:end).^(lambdaNatMen-1)*auxA./(1+THETA(indNation([2 4 7 11 16 22]))'.^2)];
ind1=find(INDACTIVE==indlambda(1)); 
ind2=[find(INDACTIVE==indNation(2)) find(INDACTIVE==indNation(4)) find(INDACTIVE==indNation(7)) find(INDACTIVE==indNation(11)) find(INDACTIVE==indNation(16)) find(INDACTIVE==indNation(22))];
se=zeros(1,length(ind2));
for i=1:length(ind2)
ind=[ind1;ind2(i)];
aux=Sigma(ind,ind);
se(i)=sqrt(D(:,i)'*aux*D(:,i));
end
Natse(1,3:end)=se;
Natse=round(ceil(Natse*1000)/1000,3);
MAT=reshape(permute(cat(3,N,Natse),[3 1 2]),16,8);
Sep=repmat([' & '],[dimNation*2 1]);
Sep1=repmat([' & ';' &('],[dimNation 1]);
Sep2=repmat(['  &  ';') & ('],[dimNation 1]);
End=repmat([' \\';')\\'],[dimNation 1]);
disp('\begin{tabular}{lcccccccc}')
disp('\hline')
disp('\hline')
disp('& \multicolumn{8}{c}{Woman}\\')
disp('&Italian& EU & New EU& Oth. Europe & Africa & Asia &S America & OECD\\')
% disp([Name1 Sep Name2 Sep num2str(N(:,1),2) Sep num2str(N(:,2),2) Sep num2str(N(:,3),2) Sep num2str(N(:,4),2) Sep num2str(N(:,5),2) Sep num2str(N(:,6),2) Sep...
%     num2str(N(:,7),2) Sep num2str(N(:,8),2) End])
disp([Name1 Sep Name2 Sep1 num2str(MAT(:,1),2) Sep2 num2str(MAT(:,2),2) Sep2 ... 
    num2str(MAT(:,3),2) Sep2 num2str(MAT(:,4),2) Sep2 num2str(MAT(:,5),2) Sep2 ... 
    num2str(MAT(:,6),2) Sep2 num2str(MAT(:,7),2) Sep2 num2str(MAT(:,8),2) End])
disp('\hline')
disp('\end{tabular}')
diary off

first=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Appendix FIGURE G1 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y=ageGrid;
X=Y';

f=figure('Name','FIGURE G1a');
plot(Y,AgeParam(2,:)','LineWidth',2)
hold on
plot(Y,AgeParam(4,:)','LineWidth',2)
hold on
plot(Y,AgeParam(6,:)','LineWidth',2)
xlabel('Age of wife')
ylabel('Effect on match surplus')
%axis([0 1 0.85 1.025])
set(gca,'fontsize',18);
l=legend('Man aged 27','Man aged 37','Man aged 50','Location','best');
set(l,'FontSize',18)
YLim=get(gca,'YLim');


f=figure('Name','FIGURE G1b');
plot(Y,AgeParam(:,2),'LineWidth',2)
hold on
plot(Y,AgeParam(:,4),'LineWidth',2)
hold on
plot(Y,AgeParam(:,6),'LineWidth',2)
xlabel('Age of husband')
ylabel('Effect on match surplus')
%axis([0 1 0.85 1.025])
set(gca,'fontsize',18);
l=legend('Woman aged 27','Woman aged 37','Woman aged 50','Location','best');
set(l,'FontSize',18)
set(gca,'YLim',YLim);
drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Appendix TABLE G1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delete('tableG1.tex')
diary('tableG1.tex')

Sep=repmat(' & ',[dimAge 1]);
End=repmat(' \\',[dimAge 1]);
prec=3;
N=round(AgeParam,prec);
disp('\begin{tabular}{llcccccc}')
disp('\hline')
disp('\hline')
disp('&& \multicolumn{6}{c}{Women}\\')
disp([' & &\multicolumn{2}{c}{Education}& & & & \\'])
disp([' & &Low & High & & & & \\'])
disp(' \cline{3-4}')
disp([' &\multicolumn{1}{r|}{Low} &' num2str(exp(A(1)),prec) '&'  num2str(exp(A(1,2)),prec)  ' &  \multicolumn{2}{c}{Housing, natives} & &  \\'])
ind=find(INDACTIVE==indA(1));aux=Sigma(ind,ind);D=A(1)*exp(A(1));aux11=sqrt(D*D*aux);
ind=find(INDACTIVE==indA(3));aux=Sigma(ind,ind);D=A(1,2)*exp(A(1,2));aux12=sqrt(D*D*aux);
disp([ '& (' num2str(round(aux11*100,2)/100,3) ')& (' num2str(round(aux12*100,2)/100,3) ') & & & & \\'])
disp(['\parbox[t]{2mm}{\multirow{3}{*}{\rotatebox[origin=c]{90}{Men}}} &\multicolumn{1}{r|}{High}&' num2str(exp(A(2,1)),prec) '&'  num2str(exp(A(2,2)),prec)  '& No & Yes  & & \\'])
ind=find(INDACTIVE==indA(3));aux=Sigma(ind,ind);D=A(2,1)*exp(A(2,1));aux21=sqrt(D*D*aux);
disp([ '& & & (' num2str(round(aux21*100,2)/100,3) ')& -'  '& &  \\'])
disp(' \cline{5-6}')
disp(['  & & &\multicolumn{1}{r|}{No}& ' num2str(1,prec) '&'  num2str(exp(lambdaH(1)),prec)  '  &\multicolumn{2}{c}{Housing, non natives} \\'])
ind=find(INDACTIVE==indlambdaH(1));aux=Sigma(ind,ind);D=lambdaH(1)*exp(lambdaH(1));aux11=sqrt(D*D*aux);
disp(['  & & & & - & ('  num2str(round(aux11*100,3)/100,prec)  ')  & & \\'])
disp(['    & &  &\multicolumn{1}{r|}{Yes}& ' num2str(exp(lambdaH(2)),prec)  '&'  num2str(exp(lambdaH(3)),prec)  '  & No &Yes  \\'])
ind=find(INDACTIVE==indlambdaH(2));aux=Sigma(ind,ind);D=lambdaH(2)*exp(lambdaH(2));aux1=sqrt(D*D*aux);
ind=find(INDACTIVE==indlambdaH(3));aux=Sigma(ind,ind);D=lambdaH(3)*exp(lambdaH(3));aux2=sqrt(D*D*aux);
disp(['    & &  & & (' num2str(round(aux1*100,3)/100,prec)  ') & ('  num2str(round(aux2*100,3)/100,prec)  '  & &  \\'])

disp(' \cline{7-8}')
disp(['    & & & & & \multicolumn{1}{r|}{No} & ' num2str(1,prec)  '&'  num2str(exp(lambdaHNI(1)),prec)  ' \\'])
ind=find(INDACTIVE==indlambdaH2(1));aux=Sigma(ind,ind);D=lambdaHNI(1)*exp(lambdaHNI(1));aux1=sqrt(D*D*aux);
disp(['    & & & & &  &  - & ('  num2str(round(aux1*100,3)/100,prec)  ') \\'])
disp(['    & & & & & \multicolumn{1}{r|}{Yes} & ' num2str(exp(lambdaHNI(2)),prec)  '&'  num2str(exp(lambdaHNI(3)),prec)  ' \\'])
ind=find(INDACTIVE==indlambdaH2(2));aux=Sigma(ind,ind);D=lambdaHNI(2)*exp(lambdaHNI(2));aux2=sqrt(D*D*aux);
ind=find(INDACTIVE==indlambdaH2(3));aux=Sigma(ind,ind);D=lambdaHNI(3)*exp(lambdaHNI(3));aux3=sqrt(D*D*aux);
disp(['    & & & & &  & (' num2str(round(aux2*100,3)/100,prec)  ') & ('  num2str(round(aux3*100,3)/100,prec)  ') \\'])
disp('\end{tabular}')
diary off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   FIGURE 10
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=figure('Name','Figure 10');
h1=scatter(Match1(:)*1000,PredictMarriage1(:)*1000,60,'MarkerEdgeColor',[0 .5 .5],'MarkerFaceColor',[0 .7 .7]);
set(gca,'fontsize',12);
hold on
h2=scatter(Match2(:)*1000,PredictMarriage2(:)*1000,60,'MarkerEdgeColor',[0 .5 .5],'MarkerFaceColor',[0.5 .1 .1]);
% title('All Marriages','FontSize',10)
xlabel('Observed Marriage Rates','fontsize',12)
ylabel('Predicted Marriage Rates','fontsize',12)
hold on
line(Match1(:)*1000,Match1(:)*1000)
set(gca,'fontsize',18);
l=legend('Before Enlargement','After Enlargement','Location','best');
drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FIGURE 12
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Culture and Access
nation=(1:8)';
auxNat=Nationality;
NoAccessM=nation>2;
NoAccessF=nation>2;
aux=logical(NoAccessM.*NoAccessF');
Beta0=ones(dimNation,dimNation);
Beta0(aux)=beta(1);
auxM=nation<=2;
auxF=nation>2;
aux=logical(auxM.*auxF');
Beta0(aux)=beta(3);
auxM=nation>2;
auxF=nation<=2;
aux=logical(auxM.*auxF');
Beta0(aux)=beta(5);
S=Nationality;
S=S.*Beta0;

NoAccessM=nation>3;
NoAccessF=nation>3;
aux=logical(NoAccessM.*NoAccessF');
Beta2=ones(dimNation,dimNation);
Beta2(aux)=beta(1);
auxM=nation<=3;
auxF=nation>3;
aux=logical(auxM.*auxF');
Beta2(aux)=beta(3);
auxM=nation>3;
auxF=nation<=3;
aux=logical(auxM.*auxF');
Beta2(aux)=beta(5);
S2=Nationality;
S2=S2.*Beta2;

SNewEU_1=S(:,3);
SNewEU_2=S2(:,3);
index=sortrows([SNewEU_2 (1:8)'],1);index=(index(:,2));


f=figure('Name','Figure 12','Position', [100, 100, 900, 1100]);
t=tiledlayout(3,1, 'TileSpacing', 'compact');

nexttile % Nationality
bh=barh([SNewEU_1(index) SNewEU_2(index)]);
yC=xline(SNewEU_1(index(end)),':k','LineWidth',3,'fontsize',18);
tC=title({'\textbf{A. Nationality Trade-offs}'});
country={'Italy','EU','EU2004-07','Other Europe','Africa','Asia','South-America','Other OECD'};
country=country(index);
set(gca,'YTickLabel',country);
a = get(gca,'YTickLabel');
set(gca,'YTickLabel',a,'FontName','Times','fontsize',16);
legC=legend(['~~Nationality' 10 '~~~+ Legal Status'],['~~Nationality' 10 '~~~~~only'],'Location','eastoutside');
legend(bh([2 1]));
legC.Interpreter = 'latex';
legC.FontSize = 20;
ylabel('Origin of Husband','FontSize',20)
xlabel('Marital Surplus','FontSize',20)
tC.Interpreter = 'latex';
tC.FontWeight='bold';
tC.FontSize=20;
tsC=subtitle('Surplus for Women from EU2004-07');
tsC.Interpreter = 'latex';
tsC.FontSize=18;

nexttile % Age
mat=AgeParam(:,2).*Nationality([1 3],3)'.*[beta(3) beta(1)];
plot(ageGrid,mat(:,1),'LineWidth',3)
hold on
plot(ageGrid,mat(:,2),'LineWidth',3)
xlabel('Age of husband','FontName','Times','fontsize',18);
ylabel({'Marriage surplus'},'FontName','Times','fontsize',18);
H=gca;
H.YLim =[0.8 2.2];
set(H,'Ygrid','on')
set(H,'Xgrid','on')
set(H,'FontSize',18)
set(H,'YMinorGrid','on')
set(H,'XMinorGrid','on')
set(H,'LineWidth',1)
 yl=yline(mat(3,2),':k','LineWidth',3);
 line([24.8;24.8],[0.8;mat(3,2)],'Color','k','LineStyle',':','LineWidth',3);
 line([39;39],[0.8;mat(3,2)],'Color','k','LineStyle',':','LineWidth',3);
 line([32;32],[0.8;mat(3,2)],'Color','k','LineStyle',':','LineWidth',3);
tA=title('\textbf{B. Age Trade-offs}');
tA.Interpreter = 'latex';
tA.FontWeight='bold';
tA.FontSize=20;
tsA=subtitle('Surplus for 27 old Women from EU2004-07');
tsA.Interpreter = 'latex';
tsA.FontSize=18;
legA = legend({['~~Surplus marrying' 10 '~~~~Italian male'],['~~Surplus marrying' 10 '~~~~EU04-07 male']});
legA.Location = 'eastoutside';
legA.Interpreter = 'latex';
legA.FontSize = 20;

nexttile % Education
mat=exp(A(:,2)).*Nationality([1 3],3)'.*[beta(3) beta(1)];
bar(mat);
xlabel('Education of husband','FontName','Times','fontsize',18);
ylabel({'Marriage surplus'},'FontName','Times','fontsize',18);
ax=gca;
% ax.YLim =[0 2];
set(gca,'XTickLabel',{'Low','High'},'FontName','Times','fontsize',18);
yl=yline(mat(2,2),':k','LineWidth',3,'fontsize',18);
tB=title('\textbf{C. Education Trade-offs}');
tB.Interpreter = 'latex';
tB.FontWeight='bold';
tB.FontSize=20;
tsB=subtitle('Surplus for High Educated Women from EU2004-07');
tsB.Interpreter = 'latex';
tsB.FontSize=18;

legB = legend({['~~Surplus marrying' 10 '~~~~Italian male'],['~~Surplus marrying' 10 '~~~~EU04-07 male']});
legB.Location = 'eastoutside';
legB.Interpreter = 'latex';
legB.FontSize = 20;
drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Solving baseline
disp('Baseline')
Glob.Policy=-1;
Period=1;
gamma0=THETA(indgamma(1));
gamma1=THETA(indgamma(2));


PredictMarriage=zeros(dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimProv*dimNum,2);

STORE=zeros(dimS,14,dimProv*dimNum);
STORE_=zeros(dimS,14,dimProv*dimNum);
GAMMA=zeros(dimProv*dimNum,1);
PROV=kron(1:dimProv,ones(1,dimNum));
NUM=kron(ones(1,dimProv),1:dimNum);
s0=zeros(2*dimS,dimProv*dimNum);
    parfor loop=1:(dimProv*dimNum)
    iprov=PROV(loop);
    num=NUM(loop);
    [Pred_Marriage,TOTAL_SURPLUS,allstore,gamma_]=solveMatch20191030(Param,Glob,iprov,Period,num);
    PredictMarriage(:,:,loop,Period)=Pred_Marriage;
    STORE(:,:,loop)=allstore.STORE;
    STORE_(:,:,loop)=allstore.STORE_;
    GAMMA(loop)=gamma_;
    s0(:,loop)=TOTAL_SURPLUS;
    end

PredictMarriage0=PredictMarriage(:,:,:,1);
PredictMarriage0=reshape(PredictMarriage0,dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimNum,dimProv);
PredictMarriage0=permute(PredictMarriage0,[1 2 4 3]);
PredictMarriage0=mean(PredictMarriage0,4); 
Singles=squeeze(STORE(:,7,:));
IdF=squeeze(STORE(:,8,:));
NatioM=squeeze(STORE(:,1,:));
NatioM_=squeeze(STORE_(:,1,:));
NatioF=squeeze(STORE(:,2,:));
NatioF_=squeeze(STORE_(:,2,:));
NatioM0=[NatioM;NatioM_];
NatioF0=[NatioF;NatioF_];

Singles_=squeeze(STORE_(:,7,:));
IdF_=squeeze(STORE_(:,8,:));

IdF0=[IdF;IdF_];
Singles0=[Singles;Singles_];

OptLove=squeeze(STORE(:,14,:));
OptLove_=squeeze(STORE_(:,14,:));
OptLove0=[OptLove;OptLove_];

% Solving counterfactual period 2 but with period 1 initial conditions
disp('Isolating the effect of policy')

Period=2;
Glob.Policy=2;
PredictMarriage=zeros(dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimProv*dimNum,2);

STORE1=zeros(dimS,14,dimProv*dimNum);
STORE1_=zeros(dimS,14,dimProv*dimNum);
PROV=kron(1:dimProv,ones(1,dimNum));
NUM=kron(ones(1,dimProv),1:dimNum);
s1=zeros(2*dimS,dimProv*dimNum);

    parfor loop=1:(dimProv*dimNum)
    iprov=PROV(loop);
    num=NUM(loop);
    [Pred_Marriage,TOTAL_SURPLUS,allstore,gamma_]=solveMatch20191030(Param,Glob,iprov,Period,num);
    PredictMarriage(:,:,loop,Period)=Pred_Marriage;
    STORE1(:,:,loop)=allstore.STORE;
    STORE1_(:,:,loop)=allstore.STORE_;
     s1(:,loop)=TOTAL_SURPLUS;     
    end

PredictMarriage1=PredictMarriage(:,:,:,1);
PredictMarriage1=reshape(PredictMarriage1,dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimNum,dimProv);
PredictMarriage1=permute(PredictMarriage1,[1 2 4 3]);
PredictMarriage1=mean(PredictMarriage1,4); 
Singles=squeeze(STORE1(:,7,:));
IdF=squeeze(STORE1(:,8,:));
Singles_=squeeze(STORE1(:,7,:));
IdF_=squeeze(STORE1_(:,8,:));
IdF1=[IdF;IdF_];
Singles1=[Singles;Singles_];
NatioM=squeeze(STORE1(:,1,:));
NatioM_=squeeze(STORE1_(:,1,:));
NatioF=squeeze(STORE1(:,2,:));
NatioF_=squeeze(STORE1_(:,2,:));
NatioM1=[NatioM;NatioM_];
NatioF1=[NatioF;NatioF_];

Glob.Policy=0;
Period=1;

auxG=[GAMMA';1-GAMMA'];
auxG=kron(auxG,ones(dimS,1))/dimS;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Surplus analysis')
SURPLUS=zeros(8,1);
OAXACA=zeros(8,4);
errlow=zeros(8,1);
errhigh=zeros(8,1);
% Italian Men married in period 1
[DSurplus,DSurplus_L,DSurplus_H,Oaxaca]=SurplusChangeMan(GAMMA,STORE,STORE_,STORE1,STORE1_,Weight,1,0);
SURPLUS(2)=DSurplus*100;
errlow(2)=DSurplus_L*100;
errhigh(2)=DSurplus_H*100;
OAXACA(2,:)=Oaxaca;

% New EU men 
[DSurplus,DSurplus_L,DSurplus_H,Oaxaca]=SurplusChangeMan(GAMMA,STORE,STORE_,STORE1,STORE1_,Weight,3,-1);
SURPLUS(4)=DSurplus*100;
errlow(4)=DSurplus_L*100;
errhigh(4)=DSurplus_H*100;
OAXACA(4,:)=Oaxaca;

% Other foreign men 
[DSurplus,DSurplus_L,DSurplus_H,Oaxaca]=SurplusChangeMan(GAMMA,STORE,STORE_,STORE1,STORE1_,Weight,-3,-1);
SURPLUS(6)=DSurplus*100;
errlow(6)=DSurplus_L*100;
errhigh(6)=DSurplus_H*100;
OAXACA(6,:)=Oaxaca;

% Italian women married in period 1
[DSurplus,DSurplus_L,DSurplus_H,Oaxaca]=SurplusChangeWoman(GAMMA,STORE,STORE_,STORE1,STORE1_,Weight,1,0);
SURPLUS(3)=DSurplus*100;
errlow(3)=DSurplus_L*100;
errhigh(3)=DSurplus_H*100;
OAXACA(3,:)=Oaxaca;

% New EU women 
[DSurplus,DSurplus_L,DSurplus_H,Oaxaca]=SurplusChangeWoman(GAMMA,STORE,STORE_,STORE1,STORE1_,Weight,3,-1);
SURPLUS(5)=DSurplus*100;
errlow(5)=DSurplus_L*100;
errhigh(5)=DSurplus_H*100;
OAXACA(5,:)=Oaxaca;

% Other foreign women 
[DSurplus,DSurplus_L,DSurplus_H,Oaxaca]=SurplusChangeWoman(GAMMA,STORE,STORE_,STORE1,STORE1_,Weight,-3,-1);
SURPLUS(7)=DSurplus*100;
errlow(7)=DSurplus_L*100;
errhigh(7)=DSurplus_H*100;
OAXACA(7,:)=Oaxaca;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FIGURE 13
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mat=[SURPLUS(2:5) 100*OAXACA(2:5,:)];

figure('Name','FIGURE 13')
bar(Mat',1)
ylabel('Percent Change in Welfare','FontSize',18)
labels = {'  Total\newlinechange in\newline welfare',...
    '  Married,\newline change in\newline   surplus',...
    '  Married,\newline change in\newline  weights',...
       ' Single\newline    to\newline married',...
       'Married\newline    to \newline single'};
a = gca;
a.XTickLabel = labels;
a = get(gca,'XTickLabel');
set(gca,'XTickLabel',a,'FontName','Times','fontsize',16)
l=legend('Native men, married','Native women, married','New EU men','New EU women','Location','best');
set(l,'FontSize',16)
drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Solving counterfactual with no EU penalty
disp('No EU penalty')
Period=1;
Glob.Policy=-1;
auxbeta0=Param.beta;
auxbeta=auxbeta0;
auxbeta(1)=1; % beta goes to 1
auxbeta(2)=0;  % no effect of shadow economy anymore
Param.beta=auxbeta;

PredictMarriage=zeros(dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimProv*dimNum,2);

STORE2=zeros(dimS,14,dimProv*dimNum);
STORE2_=zeros(dimS,14,dimProv*dimNum);
s2=zeros(2*dimS,dimProv*dimNum);
    parfor loop=1:(dimProv*dimNum)
    iprov=PROV(loop);
    num=NUM(loop);
    [Pred_Marriage,TOTAL_SURPLUS,allstore,gamma_]=solveMatch20191030(Param,Glob,iprov,Period,num);
    PredictMarriage(:,:,loop,Period)=Pred_Marriage;
    STORE2(:,:,loop)=allstore.STORE;
    STORE2_(:,:,loop)=allstore.STORE_;
    gamma(loop)=gamma_;
    end

PredictMarriage2=PredictMarriage(:,:,:,1);
PredictMarriage2=reshape(PredictMarriage2,dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimNum,dimProv);
PredictMarriage2=permute(PredictMarriage2,[1 2 4 3]);
PredictMarriage2=mean(PredictMarriage2,4); 
Singles=squeeze(STORE2(:,7,:));
IdF=squeeze(STORE2(:,8,:));
Singles_=squeeze(STORE2_(:,7,:));
IdF_=squeeze(STORE2_(:,8,:));
IdF2=[IdF;IdF_];
Singles2=[Singles;Singles_];
NatioM=squeeze(STORE2(:,1,:));
NatioM_=squeeze(STORE2_(:,1,:));
NatioF=squeeze(STORE2(:,2,:));
NatioF_=squeeze(STORE2_(:,2,:));
NatioM2=[NatioM;NatioM_];
NatioF2=[NatioF;NatioF_];

OptLove=squeeze(STORE2(:,14,:));
OptLove_=squeeze(STORE2_(:,14,:));
OptLove2=[OptLove;OptLove_];
    
    
Param.beta=auxbeta0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Solving counterfactual with surge of young African Men & Women
Glob.Policy=4;
dimAfrican=10;
sizegrid=0:5:10;
PROV=kron(1:dimProv,ones(1,dimNum));
NUM=kron(ones(1,dimProv),1:dimNum);

MM4_={zeros(dimProv*dimNum*(dimS+sizegrid(1)),14),zeros(dimProv*dimNum*(dimS+sizegrid(2)),14),zeros(dimProv*dimNum*(dimS+sizegrid(3)),14)};

for isiz=1:length(sizegrid)
disp(['Surge Young African of Size: ' num2str(dimAfrican+sizegrid(isiz))])
PolicyParam=[dimAfrican+sizegrid(isiz)];
educGrid2=kron(educGrid,ones(dimAge*dimNation,1));
educGrid3=kron(educGrid2,ones(dimHouse,1));
ageGrid2=kron(kron(ones(dimEduc,1),ageGrid),ones(dimNation,1));
ageGrid3=kron(ageGrid2,ones(dimHouse,1));
nationGrid2=kron(ones(dimAge*dimEduc,1),(1:dimNation)');
nationGrid3=kron(nationGrid2,ones(dimHouse,1));
HouseGrid3=kron(ones(dimAge*dimEduc*dimNation,1),(1:dimHouse)');


auxInd1=find(nationGrid3==5&ageGrid3==ageGrid(1)&educGrid3==10&HouseGrid3==1);
auxInd2=find(nationGrid3==5&ageGrid3==ageGrid(2)&educGrid3==10&HouseGrid3==1);
auxInd3=find(nationGrid3==5&ageGrid3==ageGrid(3)&educGrid3==10&HouseGrid3==1);
auxInd=[auxInd1;auxInd2;auxInd3];
TypeSurge=auxInd(randi(3,dimAfrican+sizegrid(isiz),dimProv*dimNum));

PredictMarriage=zeros(dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimProv*dimNum);
MM_=zeros(dimS+PolicyParam(1),8,dimProv*dimNum);

STORE4=zeros(dimS,14,dimProv*dimNum);
STORE4_=zeros(dimS+dimAfrican+sizegrid(isiz),14,dimProv*dimNum);

    parfor loop=1:(dimProv*dimNum)
    iprov=PROV(loop);
    num=NUM(loop);
    PolicyParam=[dimAfrican+sizegrid(isiz);TypeSurge(:,loop)];
    
    [Pred_Marriage,TOTAL_SURPLUS,allstore,gamma_]=solveMatch20191030(Param,Glob,iprov,Period,num,PolicyParam);
    PredictMarriage(:,:,loop,Period)=Pred_Marriage;
    STORE4(:,:,loop)=allstore.STORE;
    STORE4_(:,:,loop)=allstore.STORE_;
    gamma(loop)=gamma_;
    end

PredictMarriage4=mean(PredictMarriage,3);
MM4_{isiz}=reshape(permute(STORE4_,[1 3 2]),(dimS+PolicyParam(1))*dimProv*dimNum,14);

end

Glob.Policy=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Graph

Married=zeros(3,2,length(sizegrid)+1);
AgeGap=zeros(3,2,length(sizegrid)+1);
%                Baseline %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
auxM=reshape(permute(STORE_,[1 3 2]),dimS*dimProv*dimNum,14);
% African men, not single, married with Italian women
ind2=auxM(:,1)==5&(auxM(:,3)==27|auxM(:,3)==32); 
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)==1&(auxM(:,3)<=32); 
Married(1,1,1)=sum(ind)/sum(ind2);
% African men, not single, married with African women
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)==5&(auxM(:,3)<=32); 
Married(2,1,1)=sum(ind)/sum(ind2);
% African men, not single, married with any women
ind=auxM(:,7)==0&auxM(:,1)==5&(auxM(:,3)<=32); 
Married(3,1,1)=sum(ind)/sum(ind2);
% African women, not single, married with Italian men
ind2=auxM(:,2)==5&(auxM(:,4)==27|auxM(:,4)==32); 
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)==1&(auxM(:,4)<=32); 
Married(1,2,1)=sum(ind)/sum(ind2);
% African men, not single, married with African women
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)==5&(auxM(:,4)<=32); 
Married(2,2,1)=sum(ind)/sum(ind2);
% African men, not single, married with any women
ind=auxM(:,7)==0&auxM(:,2)==5&(auxM(:,4)<=32); 
Married(3,2,1)=sum(ind)/sum(ind2);

% African men, not single, married with Italian women
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)==1&(auxM(:,3)<=32); 
AgeGap(1,1,1)=mean(auxM(ind,3)-auxM(ind,4)); 
% African men, not single, married with African women
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)==5&(auxM(:,3)<=32); 
AgeGap(2,1,1)=mean(auxM(ind,3)-auxM(ind,4)); 
% African men, not single, married with other women
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)~=1&auxM(:,2)~=5&(auxM(:,3)<=32); 
AgeGap(3,1,1)=mean(auxM(ind,3)-auxM(ind,4)); 

% African women, not single, married with Italian men
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)==1&(auxM(:,4)<=32); 
AgeGap(1,2,1)=mean(auxM(ind,3)-auxM(ind,4)); % Not single, married with Italian
% African women, not single, married with African men
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)==5&(auxM(:,4)<=32); 
AgeGap(2,2,1)=mean(auxM(ind,3)-auxM(ind,4)); 
% African women, not single, married with other men
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)~=1&auxM(:,1)~=5&(auxM(:,4)<=32); 
AgeGap(3,2,1)=mean(auxM(ind,3)-auxM(ind,4)); 


for isiz=1:length(sizegrid)
    auxM=MM4_{isiz};

% African men, not single, married with Italian women
ind2=auxM(:,1)==5&(auxM(:,3)==27|auxM(:,3)==32); 
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)==1&(auxM(:,3)<=32); 
Married(1,1,isiz+1)=sum(ind)/sum(ind2);
% African men, not single, married with African women
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)==5&(auxM(:,3)<=32); 
Married(2,1,isiz+1)=sum(ind)/sum(ind2);
% African men, not single, married with any women
ind=auxM(:,7)==0&auxM(:,1)==5&(auxM(:,3)<=32); 
Married(3,1,isiz+1)=sum(ind)/sum(ind2);

% African women, not single, married with Italian men
ind2=auxM(:,2)==5&(auxM(:,4)==27|auxM(:,4)==32); 
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)==1&(auxM(:,4)<=32); 
Married(1,2,isiz+1)=sum(ind)/sum(ind2);
% African men, not single, married with African women
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)==5&(auxM(:,4)<=32); 
Married(2,2,isiz+1)=sum(ind)/sum(ind2);
% African men, not single, married with any women
ind=auxM(:,7)==0&auxM(:,2)==5&(auxM(:,4)<=32); 
Married(3,2,isiz+1)=sum(ind)/sum(ind2);

% African men, not single, married with Italian women
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)==1&(auxM(:,3)<=32); 
AgeGap(1,1,isiz+1)=mean(auxM(ind,3)-auxM(ind,4)); 
% African men, not single, married with African women
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)==5&(auxM(:,3)<=32); 
AgeGap(2,1,isiz+1)=mean(auxM(ind,3)-auxM(ind,4)); 
% African men, not single, married with other women
ind=auxM(:,7)==0&auxM(:,1)==5&auxM(:,2)~=1&auxM(:,2)~=5&(auxM(:,3)<=32); 
AgeGap(3,1,isiz+1)=mean(auxM(ind,3)-auxM(ind,4)); 

% African women, not single, married with Italian men
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)==1&(auxM(:,4)<=32); 
AgeGap(1,2,isiz+1)=mean(auxM(ind,3)-auxM(ind,4)); % Not single, married with Italian
% African women, not single, married with African men
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)==5&(auxM(:,4)<=32); 
AgeGap(2,2,isiz+1)=mean(auxM(ind,3)-auxM(ind,4)); 
% African women, not single, married with other men
ind=auxM(:,7)==0&auxM(:,2)==5&auxM(:,1)~=1&auxM(:,1)~=5&(auxM(:,4)<=32); 
AgeGap(3,2,isiz+1)=mean(auxM(ind,3)-auxM(ind,4)); 

end
Married=permute(Married,[1 3 2]);
AgeGap=permute(AgeGap,[1 3 2]);
Married=Married(:,[1 2 4],:);
AgeGap=AgeGap(:,[1 2 4],:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   FIGURE 15
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('Name','FIGURE 15','Position', [100, 100, 1000, 1000])
t=tiledlayout(2,2, 'TileSpacing', 'compact');

nexttile % Men
b1=bar(Married(:,:,1),0.5);
ylabel('Percent married')
labels = {'     With\newlineItalian spouse','     With\newlineAfrican spouse','  Any\newline spouse'};
a1 = gca;
a1.XTickLabel = labels;
a11 = get(gca,'XTickLabel');
set(gca,'XTickLabel',a11,'FontName','Times','fontsize',16)
t1=title('\textbf{Married, African Men}');
t1.Interpreter = 'latex';
ts1=subtitle('Age $\le$ 32');
ts1.Interpreter = 'latex';

l=legend('Baseline','7% increase','15% increase','Location','eastoutside');
set(l,'FontSize',18)
title(l,'Size of African Community')
l.Layout.Tile = 'South';
l.Orientation='Horizontal';
nexttile % Women
b2=bar(Married(:,:,2),0.5);
ylabel('Percent married')
labels = {'     With\newlineItalian spouse','     With\newlineAfrican spouse','  Any\newline spouse'};
a2 = gca;
a2.XTickLabel = labels;
a22 = get(gca,'XTickLabel');
set(gca,'XTickLabel',a22,'FontName','Times','fontsize',16)
t2=title('\textbf{Married, African Women}');
t2.Interpreter = 'latex';
ts2=subtitle('Age $\le$ 32');
ts2.Interpreter = 'latex';

a1.YLim=a2.YLim;

nexttile % Men
b1=bar(AgeGap(:,:,1),0.5);
ylabel('Average Age Gap (Years)')
labels = {'     With\newlineItalian spouse','     With\newlineAfrican spouse','  Other\newline spouse'};
a3 = gca;
a3.XTickLabel = labels;
a11 = get(gca,'XTickLabel');
set(gca,'XTickLabel',a11,'FontName','Times','fontsize',16)
t3=title('\textbf{Age Gap, African Men}');
t3.Interpreter = 'latex';
ts3=subtitle('Age $\le$ 32');
ts3.Interpreter = 'latex';

nexttile(4) % Men
b1=bar(AgeGap(:,:,2),0.5);
ylabel('Average Age Gap (Years)')
labels = {'     With\newlineItalian spouse','     With\newlineAfrican spouse','  Other\newline spouse'};
a4 = gca;
a4.XTickLabel = labels;
a11 = get(gca,'XTickLabel');
set(gca,'XTickLabel',a11,'FontName','Times','fontsize',16)
t4=title('\textbf{Age Gap, African Women}');
t4.Interpreter = 'latex';
ts4=subtitle('Age $\le$ 32');
ts4.Interpreter = 'latex';

a3.YLim=[min(AgeGap(:)) ceil(max(AgeGap(:))) ];
a4.YLim=[min(AgeGap(:)) ceil(max(AgeGap(:)))];
drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Solving counterfactual with surge balanced in terms of nationality
Glob.Policy=5;
dimNewMigrant=20;
PolicyParam=[dimNewMigrant];
MM5_=zeros(dimProv*dimNum*(dimS+PolicyParam(1)),14,2);
Origin=5;

for ageloop=2:2
    disp(['Surge Young Migrants of age : ' num2str(ageGrid(ageloop))])
    PolicyParam=[dimNewMigrant;ageGrid(ageloop);Origin];
    Glob.PolicyParam=PolicyParam;
    PredictMarriage=zeros(dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimProv*dimNum);
    MM_=zeros(dimS+PolicyParam(1),8,dimProv*dimNum);
    STORE5=zeros(dimS,14,dimProv*dimNum);
    STORE5_=zeros(dimS+dimNewMigrant,14,dimProv*dimNum);

    parfor loop=1:(dimProv*dimNum)
    iprov=PROV(loop);
    num=NUM(loop);
    [Pred_Marriage,TOTAL_SURPLUS,allstore,gamma_]=solveMatch20191030(Param,Glob,iprov,Period,num);
    PredictMarriage(:,:,loop,Period)=Pred_Marriage;
    STORE5(:,:,loop)=allstore.STORE;
    STORE5_(:,:,loop)=allstore.STORE_;
    gamma(loop)=gamma_;
    end

PredictMarriage4=mean(PredictMarriage,3);
MM5_(:,:,ageloop)=reshape(permute(STORE5_,[1 3 2]),(dimS+PolicyParam(1))*dimProv*dimNum,14);

end

Glob.Policy=0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Displaying results

% IndForeign=not((repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==1)|(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])==1));
IndForeign=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==repmat(nationGrid,[1 dimAge*dimEduc*dimNation])')&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])~=1);
IndForeign=repmat(IndForeign,[1 1 dimProv]);
Ind1_1=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==1)&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])==1);
Ind1_1=repmat(Ind1_1,[1 1 dimProv]);
HeterM=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==1)&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])>=2);
HeterF=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])>=2)&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])==1);
Heter=HeterM==1|HeterF==1;
Heter=repmat(Heter,[1 1 dimProv]);

aux=mean(PredictMarriage0,3);
total_0=sum(aux(:));
aux=mean(PredictMarriage2,3);
total_2=sum(aux(:));
total=(total_2-total_0)/total_0;

aux=mean(PredictMarriage0(Ind1_1),3);
q11_0=sum(aux(:));
aux=mean(PredictMarriage2(Ind1_1),3);
q11_2=sum(aux(:));
Ita_Ita=(q11_2-q11_0)/q11_0;

aux=mean(PredictMarriage0(Heter),3);
qHeter_0=sum(aux(:));
aux=mean(PredictMarriage2(Heter),3);
qHeter_2=sum(aux(:));
Ita_Other=(qHeter_2-qHeter_0)/qHeter_0;

aux=mean(PredictMarriage0(IndForeign),3);
Foreign0=sum(aux(:));
aux=mean(PredictMarriage2(IndForeign),3);
Foreign2=sum(aux(:));
Foreigners=(Foreign2-Foreign0)/Foreign0;

ufertility=randn(size(OptLove0));
fertility0=double((gamma0+gamma1*OptLove0+ufertility)>0);
fertility0(logical(Singles0))=NaN;


fertility2=double((gamma0+gamma1*OptLove2+ufertility)>0);
fertility2(logical(Singles2))=NaN;

Ind0=1-Singles0;Ind0(Singles0==1)=NaN;
fertility0=normcdf(gamma0+gamma1*OptLove0.*Ind0);

Ind2=1-Singles2;Ind2(Singles2==1)=NaN;
fertility2=normcdf(gamma0+gamma1*OptLove2.*Ind2);


% Fertility Homogamous Natives
auxF0_1=mean(nanmean(fertility0(1:dimS,:))); % in the first pool
Ind=NatioM0(dimS+1:end,:)==1&NatioF0(dimS+1:end,:)==1;
auxF0_2=fertility0(dimS+1:end,:);
auxF0_2=nanmean(auxF0_2(Ind)); % in the second pool
auxF0=mean(GAMMA)*auxF0_1+(1-mean(GAMMA))*auxF0_2;

auxF2_1=mean(nanmean(fertility2(1:dimS,:))); % in the first pool
Ind=NatioM2(dimS+1:end,:)==1&NatioF2(dimS+1:end,:)==1;
auxF2_2=fertility2(dimS+1:end,:);
auxF2_2=nanmean(auxF2_2(Ind)); % in the second pool
auxF2=mean(GAMMA)*auxF2_1+(1-mean(GAMMA))*auxF2_2;
DFert_Ita_Ita=auxF2/auxF0-1;

% Fertility Heterogamous
Ind=(NatioM0(dimS+1:end,:)~=NatioF0(dimS+1:end,:))&((NatioM0(dimS+1:end,:)==1)|(NatioF0(dimS+1:end,:)==1));
auxF0=fertility0(dimS+1:end,:);
auxF0=nanmean(auxF0(Ind)); % in the second pool
% Ind=(NatioM2(dimS+1:end,:)~=NatioF2(dimS+1:end,:));
Ind=(NatioM2(dimS+1:end,:)~=NatioF2(dimS+1:end,:))&((NatioM2(dimS+1:end,:)==1)|(NatioF2(dimS+1:end,:)==1));
auxF2=fertility2(dimS+1:end,:);
auxF2=nanmean(auxF2(Ind)); % in the second pool
DFert_Heter=auxF2/auxF0-1;

% Fertility Homogamous Foreigners
Ind=(NatioM0(dimS+1:end,:)==NatioF0(dimS+1:end,:))&NatioM0(dimS+1:end,:)~=1;
auxF0=fertility0(dimS+1:end,:);
auxF0=nanmean(auxF0(Ind)); % in the second pool
Ind=(NatioM2(dimS+1:end,:)==NatioF2(dimS+1:end,:))&NatioM2(dimS+1:end,:)~=1;
auxF2=fertility2(dimS+1:end,:);
auxF2=nanmean(auxF2(Ind)); % in the second pool
DFert_Hom_For=auxF2/auxF0-1;

% All marriages 
auxF0_1=mean(nanmean(fertility0(1:dimS,:))); % in the first pool
auxF0_2=mean(nanmean(fertility0(dimS+1:end,:))); % in the second pool
auxF0=mean(GAMMA)*auxF0_1+(1-mean(GAMMA))*auxF0_2;

auxF2_1=mean(nanmean(fertility2(1:dimS,:))); % in the first pool
auxF2_2=mean(nanmean(fertility2(dimS+1:end,:))); % in the second pool
auxF2=mean(GAMMA)*auxF2_1+(1-mean(GAMMA))*auxF2_2;
DFert_All=auxF2/auxF0-1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  FIGURE 14
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=figure('Name','Figure 14','Position', [100, 100, 1900, 1000]);
t=tiledlayout(1,2, 'TileSpacing', 'compact');

nexttile % Marriage
qMarriage=[Foreigners;Ita_Other;Ita_Ita; total]*100;
barh(qMarriage,0.5)
xlabel('Percentage change in marriage rates')
labels = {'   All\newlinemarriages','Homogamous\newline  marriages\newline   (natives)','Heterogamous\newline   marriages','Homogamous\newline   marriages\newline (foreigners)'};
labels = {'   All marriages','Homogamous marriages\newline  (natives)','Heterogamous marriages\newline  (with natives)','Homogamous marriages\newline (foreigners)'};
labels = {'Homogamous marriages\newline (foreigners)','Heterogamous marriages\newline  (with natives)','Homogamous marriages\newline  (natives)','   All marriages'};
a = gca;
a.YTickLabel = labels;
a = get(gca,'YTickLabel');
set(gca,'YTickLabel',a,'FontName','Times','fontsize',30)
tC1=title({'\textbf{Marriages}'});
tC1.Interpreter = 'latex';
tC1.FontWeight='bold';
tC1.FontSize=30;

nexttile;
qFert=[DFert_Hom_For;DFert_Heter;DFert_Ita_Ita; DFert_All]*100;
barh(qFert,0.5)
xlabel('Percentage change in fertility rates')
labels = {' ',' ',' ',' '};
a = gca;
a.YTickLabel = labels;
a = get(gca,'YTickLabel');
set(gca,'YTickLabel',a,'FontName','Times','fontsize',30)
tC2=title({'\textbf{Fertility}'});
tC2.Interpreter = 'latex';
tC2.FontWeight='bold';
tC2.FontSize=30;
drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FIGURE 11
dimNum=1000;
dimProv=21;


%% Solving baseline
disp('Baseline for FIGURE 11')
tic
Glob.Policy=-1;
PolicyParam=0;
Glob.PolicyParam=PolicyParam;

Period=1;

PredictMarriage=zeros(dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimProv*dimNum);

STORE=zeros(dimS,14,dimProv*dimNum);
STORE_=zeros(dimS,14,dimProv*dimNum);
GAMMA=zeros(dimProv*dimNum,1);
PROV=kron(1:dimProv,ones(1,dimNum));
NUM=kron(ones(1,dimProv),1:dimNum);
gamma=zeros(dimProv*dimNum,1);

    parfor loop=1:(dimProv*dimNum)
    iprov=PROV(loop);
    num=NUM(loop);
    [Pred_Marriage,TOTAL_SURPLUS,allstore0,gamma_,Surplus0,Surplus0_]=solveMatch20191030(Param,Glob,iprov,Period,num);
    PredictMarriage(:,:,loop,Period)=Pred_Marriage;
    STORE(:,:,loop)=allstore0.STORE;
    STORE_(:,:,loop)=allstore0.STORE_;
    gamma(loop)=gamma_;
    end

PredictMarriage0=reshape(PredictMarriage,dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimNum,dimProv);
PredictMarriage0=permute(PredictMarriage0,[1 2 4 3]);
PredictMarriage0=mean(PredictMarriage0,4); 

Singles0=squeeze(STORE(:,7,:));
Singles0=Singles0(:);
typeM=squeeze(ceil(STORE(:,12,:)/2));
typeM=typeM(:);
typeF=squeeze(ceil(STORE(:,13,:)/2));
typeF=typeF(:);

SingleTypeM=typeM(logical(Singles0));
SingleTypeF=typeF(logical(Singles0));
aux=tabulate(SingleTypeM);
TABM0=zeros(dimAge*dimEduc*dimNation,1);
TABM0_=zeros(dimAge*dimEduc*dimNation,1);
TABM0(1:size(aux,1),1)=aux(:,2);
aux=tabulate(SingleTypeF);
TABF0=zeros(dimAge*dimEduc*dimNation,1);
TABF0_=zeros(dimAge*dimEduc*dimNation,1);
TABF0(1:size(aux,1),1)=aux(:,2);


Singles0_=squeeze(STORE_(:,7,:));
Singles0_=[Singles0_;nan(1,dimNum*dimProv)];
Singles0_=reshape(Singles0_,(dimS+1)*dimNum*dimProv,1);

typeM0_=squeeze(ceil(STORE_(:,12,:)/2));
typeM0_=typeM0_(:);
typeF0_=squeeze(ceil(STORE_(:,13,:)/2));
typeF0_=typeF0_(:);
OptF0_=squeeze(STORE_(:,8,:));
OptF0_=[OptF0_;nan(1,dimNum*dimProv)];
OptF0_=reshape(OptF0_,(dimS+1)*dimNum*dimProv,1);

SingleTypeM_=typeM0_(logical(Singles0_(not(isnan(Singles0_)))));
SingleTypeF_=typeF0_(logical(Singles0_(not(isnan(Singles0_)))));
aux=tabulate(SingleTypeM_);
TABM0_(1:size(aux,1),1)=aux(:,2);
aux=tabulate(SingleTypeF_);
TABF0_(1:size(aux,1),1)=aux(:,2);

auxgamma=mean(gamma);
BaseM=TABM0*auxgamma+(1-auxgamma)*TABM0_;
BaseF=TABF0*auxgamma+(1-auxgamma)*TABF0_;
BaseM=repmat(BaseM,[1 dimAge*dimEduc*dimNation]);
BaseF=repmat(BaseF,[1 dimAge*dimEduc*dimNation]);
toc
Glob.Policy=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Period=1;
Glob.Policy=6;
DistMen=Glob.Singles_Men_Before;
aux=sum(DistMen,2);
WeightElast=aux/max(aux);
WeightElast=repmat(WeightElast,[1 dimAge*dimEduc*dimNation]);

educGrid2=kron(educGrid,ones(dimAge*dimNation,1));
educGrid3=kron(educGrid2,ones(dimHouse,1));
ageGrid2=kron(kron(ones(dimEduc,1),ageGrid),ones(dimNation,1));
ageGrid3=kron(ageGrid2,ones(dimHouse,1));
nationGrid2=kron(ones(dimAge*dimEduc,1),(1:dimNation)');
nationGrid3=kron(nationGrid2,ones(dimHouse,1));
HouseGrid3=kron(ones(dimAge*dimEduc*dimNation,1),(1:dimHouse)');

STORE4=zeros(dimS+1,14,dimProv*dimNum,96);
STORE4_=zeros(dimS+1,14,dimProv*dimNum,96);
PredictMarriage=zeros(dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimProv*dimNum,96);

for bigloop=1:length(nationGrid)
disp(['Men of age : ' num2str(ageGrid2(bigloop)) ', educ : ' num2str(educGrid2(bigloop)) ', origin : ' num2str(nationGrid(bigloop))  ])
PolicyParam=[bigloop*2];  % *2 as we pick individuals with a house
Glob.PolicyParam=PolicyParam;

    parfor loop=1:(dimProv*dimNum)
    iprov=PROV(loop);
    num=NUM(loop);
    [Pred_Marriage,TOTAL_SURPLUS,allstore,gamma_,Surplus,Surplus_]=solveMatch20191030(Param,Glob,iprov,Period,num);
        PredictMarriage(:,:,loop,bigloop)=Pred_Marriage;
    STORE4(:,:,loop,bigloop)=allstore.STORE;    
    STORE4_(:,:,loop,bigloop)=allstore.STORE_;    
    end

end

STORE44=permute(STORE4,[4 2 1 3]); % dPol(=96) 10 (dimS) dProv*dNum
STORE44=reshape(STORE44,96,14,(dimS+1)*dimNum*dimProv);
STORE44=permute(STORE44,[3 2 1]); % dSim 10 dPol
typeM=squeeze(ceil(STORE44(:,12,:)/2));
typeF=squeeze(ceil(STORE44(:,13,:)/2));

Singles=STORE44(:,7,:);
Singles=permute(Singles,[1 3 2]);

STORE44_=permute(STORE4_,[4 2 1 3]); % dPol(=96) 10 (dimS) dProv*dNum
STORE44_=reshape(STORE44_,96,14,(dimS+1)*dimNum*dimProv);
STORE44_=permute(STORE44_,[3 2 1]); % dSim 10 dPol
typeM_=squeeze(ceil(STORE44_(:,12,:)/2));
typeF_=squeeze(ceil(STORE44_(:,13,:)/2));
OptF_=squeeze(STORE44_(:,8,:));
Singles_=STORE44_(:,7,:);
Singles_=permute(Singles_,[1 3 2]);

TABM=zeros(96,96); % Number of singles men by type across all markets
TABF=zeros(96,96);
TABM_=zeros(96,96);
TABF_=zeros(96,96);

for bigloop=1:(dimAge*dimEduc*dimNation)
    Ikeep=not(isnan(Singles(:,bigloop)));
    auxSingle=logical(Singles(Ikeep,bigloop));
    auxTypeM=typeM(Ikeep,bigloop);
    SingleTypeM=auxTypeM(auxSingle,1);
    auxTypeF=typeF(Ikeep,bigloop);
    SingleTypeF=auxTypeF(auxSingle,1);
    aux=tabulate(SingleTypeM);
    TABM(1:length(aux),bigloop)=aux(:,2);
    aux=tabulate(SingleTypeF);
    TABF(1:length(aux),bigloop)=aux(:,2);

    Ikeep=logical(Singles(:,bigloop));
    SingleTypeM=typeM(Ikeep,bigloop);
    SingleTypeF=typeF(Ikeep,bigloop);
    aux=tabulate(SingleTypeM);
    TABM(1:size(aux,1),bigloop)=aux(:,2);
    aux=tabulate(SingleTypeF);
    TABF(1:size(aux,1),bigloop)=aux(:,2);


    Ikeep=logical(Singles_(:,bigloop));
    SingleTypeM_=typeM_(Ikeep,bigloop);
    SingleTypeF_=typeF_(Ikeep,bigloop);
    aux=tabulate(SingleTypeM_);
    TABM_(1:size(aux,1),bigloop)=aux(:,2);
    aux=tabulate(SingleTypeF_);
    TABF_(1:size(aux,1),bigloop)=aux(:,2);
    TABF(1,bigloop)=TABF(1,bigloop)-dimNum*dimProv; % This is the fictious woman introduced to make the market square
    TABF_(1,bigloop)=TABF_(1,bigloop)-dimNum*dimProv; % This is the fictious woman introduced to make the market square
    TABM(1,bigloop)=TABM(1,bigloop)-dimNum*dimProv*(nationGrid(bigloop)~=1); % This is the fictious man introduced to make the market square
end

Glob.Policy=0;
Glob.PolicyParam=0;


DistMen=Glob.Singles_Men_Before;
CounterM=TABM*auxgamma+(1-auxgamma)*TABM_;
CounterF=TABF*auxgamma+(1-auxgamma)*TABF_;
CounterM=WeightElast.*CounterM+(1-WeightElast).*BaseM;
CounterF=WeightElast.*CounterF+(1-WeightElast).*BaseF;

indExtract=tril(true(dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation),-1);
RLK_MM=CounterM./BaseM-1;
YM=RLK_MM(indExtract);
aux=RLK_MM';
XM=aux(indExtract);

RLK_FM=CounterF./BaseF-1;
YF=RLK_FM(indExtract);
aux=RLK_FM';
XF=aux(indExtract);


MATXY=sortrows([XM YM],1);
ind=(sum(isnan(MATXY),2)==0)&(sum(isinf(MATXY),2)==0);
MATXY=MATXY(ind,:);
X=[ones(length(MATXY),1) MATXY(:,1)];Y=MATXY(:,2);
b=inv(X'*X)*X'*Y;
predY=X*b;
predY=csaps(X(:,2),Y,0.15,X(:,2));

f=figure('Name','FIGURE 11a');
scatter(MATXY(:,1),MATXY(:,2))
hold on
plot(MATXY(MATXY(:,1)>0,1),MATXY(MATXY(:,1)>0,1),'LineWidth',2)  % Choo & Siow
xlabel('R_{xx\prime}')
ylabel('R_{x\prime x}')
set(gca,'fontsize',18);
l=legend('Model','Choo & Siow (45^o degree line)','Location','best');
set(l,'FontSize',18)
YLim=get(gca,'YLim');
xlim([0 Inf])
ylim([0 Inf])


MATXY=sortrows([XF YF],1);
ind=(sum(isnan(MATXY),2)==0)&(sum(isinf(MATXY),2)==0);
MATXY=MATXY(ind,:);
X=[ones(length(MATXY),1) MATXY(:,1)];Y=MATXY(:,2);
b=inv(X'*X)*X'*Y;
predY=X*b;
predY=csaps(X(:,2),Y,0.85,X(:,2));

f=figure('Name','FIGURE 11b');
scatter(MATXY(:,1),MATXY(:,2))
hold on
plot(MATXY(MATXY(:,1)<0,1),MATXY(MATXY(:,1)<0,1),'LineWidth',2)  % Choo & Siow
xlabel('R_{xx\prime}')
ylabel('R_{x\prime x}')
set(gca,'fontsize',18);
l=legend('Model','Choo & Siow (45^o degree line)','Location','best');
set(l,'FontSize',18)
xlim([-Inf 0])
ylim([-Inf 0])
YLim=get(gca,'YLim');

end

if stderr==1
    disp('Computing Standard Errors')
    dimMomentsM=2*dimAge*dimEduc*dimNation*dimAge*dimEduc*dimNation*dimProv;
    dimMoments=dimMomentsM+7*3+4;
INDACTIVE=[indOutside(1:22) indA(1:3) indAgeParam inddelta indsigma indNation indbeta indlambdaH2 indlambda([1 4 5]) indlambdaH(1:3) indgamma(1:2)];


    indactive=INDACTIVE;
    Sstore=zeros(dimMoments,length(indactive));
    not_done=(sum(Sstore)==0)|isnan(sum(Sstore));
    prec=0.00001;
   
    iter=1;
    while sum(not_done)~=0
        
       disp(['Iter = ' num2str(iter) '      # of parameters =  ' num2str(sum(not_done==1))  '   Precision = ' num2str(prec)])
       indactive=INDACTIVE(not_done);
       S=gradp_forward(@estimMarriage20191030,THETA(indactive),prec,'forward3');
       Sstore(:,not_done)=S;
       disp(num2str(sum(Sstore)',5))
       not_done=(sum(abs(Sstore))==0)|isnan(sum(Sstore)); 
       prec=prec*10;
       iter=iter+1;
       if prec==100;
           break
       end
    end
    H=Sstore;
    auxD=DivorceMoments(:,2,:);
    denom=[ones(dimMomentsM,1)/10;auxD(:)];
    denom=[SD_Match1(:);SD_Match2(:);auxD(:);ObsFertility(2,:)'];
    Omega=speye(dimMoments);
    Omega(logical(Omega))=denom;
    

    Sigma=pinv(H'/Omega*H);
    num2str(sqrt(diag(Sigma)))
    setheta=nan(length(THETA),1);
    setheta(INDACTIVE)=sqrt(diag(Sigma));
      save Sstore2 Sstore Sigma INDACTIVE setheta
    
 
    
end
